
lagrange := proc(a::list, b::list)
local i, qi, n, P, A, Ai;

    P := 0:
    n := nops(a):

    A := expand(mul(X - a[k], k = 1..n)):

    for i to n do
        Ai := simplify(A / (X - a[i])):

        qi := subs(X = a[i], Ai):

        P := P + b[i] * Ai / qi:
    end do:

    expand(P);
end proc;
